Hydrodynamic Stability of Cosmological Quark-Hadron Phase Transitions 



P. Chris Fragile and Peter Anninos 
University of California, Lawrence Livermore National Laboratory, Livermore, CA 94550 

Abstract 

Recent results from linear perturbation theory suggest that first-order cosmological quark-hadron phase transitions occurring 
as deflagrations may be "borderline" unstable, and those occurring as detonations may give rise to growing modes behind the 
interface boundary. However, since nonlinear effects can play important roles in the development of perturbations, unstable 
behavior cannot be asserted entirely by linear analysis, and the uncertainty of these recent studies is compounded further 
by nonlinearities in the hydrodynamics and self-interaction fields. In this paper we investigate the growth of perturbations 
and the stability of quark-hadron phase transitions in the early Universe by solving numerically the fully nonlinear relativistic 
hydrodynamics equations coupled to a scalar field with a quartic self-interaction potential regulating the transitions. We 
consider single, perturbed, phase transitions propagating either by detonation or deflagration, as well as multiple phase and 
shock front interactions in 1+2 dimensional spacetimes. 
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I. INTRODUCTION 

First-order phase transitions ocurring at either the electroweak or QCD symmetry breaking epochs in the early 
Universe, as predicted by the standard model of cosmology, can have important consequences for the history of our 
Universe. In particular, the formation of co-existing bubbles and droplets of different phases during the QCD transition 
affects the production and distribution of hadrons, and may lead to significant baryon number density fluctuations. 
Assuming these fluctuations survive to the epoch of Big Bang nucleosynthesis, they can strongly influence the predicted 
light element abundances 0. This, in turn, may alter our current understanding and interpretation of homogeneous 
Big Bang nucleosynthesis, as well as the dynamical and chemical evolution of matter structures, including the origins 
of primordial magnetic fields and structure perturbations that seed the subsequent formation of stellar, galactic, and 
cluster scale systems. 

Baryon density fluctuations evolve hydrodynamically by the competing effects of local phase mixing and phase 
separation that may arise during the transition period. A potential mixing mechanism that we consider here is the 
hydrodynamic instability of the interface surfaces separating regions of different phase. Unstable modes may distort the 
phase boundary and induce mixing and diffusive homogenization through hydrodynamic turbulence. Whether unstable 
modes can exist in cosmological phase fronts has been discussed recently by several authors with conflicting conclusions 
based on first-order perturbation theory. For example, Kamionkowski and Freese [2f suggest that subsonic deflagration 
fronts in electroweak transitions may be accelerated until they become supersonic and proceed as detonations through 
an effective increase in the front velocity due to surface distortion effects as the transition becomes turbulent. Link 
indicated that the phase boundary in slow deflagrations from quark-hadron transitions may be unstable to long 
wavelength perturbations, but may be stabilized by surface tension at shorter wavelengths. In contrast, Huet et al. |4| 
find that for electroweak transitions, the shape of the phase boundary is stable under perturbations, and quark-hadron 
transitions are at the "borderline" between stable and unstable. They also argue that unstable modes are not possible 
in detonation fronts. However, Abney [f| suggests that this is the case for the quark phase ahead of the supersonic 
interface boundary, but that the cold phase hadron regions behind the bubbles might be unstable, at least for the 
Chapman- Jouget type detonations investigated. 

These results are certainly not conclusive, and in some cases are in apparent contradiction due to the level of 
approximation and coupling assumed between thermal and dynamical variables. Also, the full consequences or even 
existence of instabilities cannot be determined entirely by linear analysis. Nonlinear effects, including higher-order 
coupling between hydrodynamic, microphysical, scalar field, and interaction potentials, as well as surface tension, 
entropy flow, and baryon transport, may all play important roles in the stabilization (or de-stabilization) of phase 
boundaries, and remain to be investigated in detail. Hence, we undertake this study to more fully investigate the 
hydrodynamic stability of cosmological phase transitions ocurring during the QCD epoch. We take a numerical 
approach and thus generalize previous results by solving the multi-dimensional relativistic hydrodynamics equations, 
presented in EflTl coupled together with a model scalar wave equation and an interaction potential regulating the phase 
transitions. Using results from perturbation theory to define initial data (as discussed in 31110 . we solve numerically 
for the evolutions of single distorted phase fronts as well as interactions and collisions of multiple front systems 
propagating as either deflagrations or detonations. Our numerical results are presented in H1VI and summarized in 



II. BASIC EQUATIONS 

The equations formulated with internal fluid energy are derived from the 4- velocity normalization u^u^ = — 1, 
the parallel component of the stress-energy conservation equation m„V = for internal energy, the transverse 
component of the stress-energy equation (g Q „ + u Q w„)V M T Ml ' = for momentum, and an equation of state for the fluid 
pressure and temperature, using traditional high energy units in which c = h = kg = 1. We consider the following 
special relativistic stress-energy tensor 



which includes both fluid and scalar field (</>) contributions, ph is the relativistic enthalpy, p r is the radiation pressure, 
is the contravariant 4-velocity, and g M „ = 7^ t „ is the flat space 4-metric. 
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FIG. 1: Interaction potential V(<f>, T) (equation i'2.'2h in the text) as a function of scalar field at four different temperatures 
T/T c — (2, 1, 0.9, 0.5). The potential parameters used here are: surface tension a = 0.1 T c , correlation length £ c — 1.0/T C , 
latent heat L — 2.0 , and critical temperature T c = 1. 



where T is the fluid temperature, and A, a, /3, and To are constant parameters specified indirectly through their 
relationships to the critical transition temperature T c , surface tension a, correlation length £ c , and latent heat L 0,0 



To , a \ _ f T c [l + da/iLQ]- 1 ' 2 , (2a^T c 2 /Z)- 1 ' 2 
P, A J " \ (L + 6a/£ c ) (QalcT?)- 1 , (3a^)-i 



(2.3) 



The potential 11'. 21 is plotted at four different temperatures T/T c = (2, 1, 0.9, 0.5) as a function of scalar field <f> 
in Figure ^ At high temperatures (T > T c ) the potential exhibits a single minimum (0 — V{4>, T) — 0), where 
the cosmological fluid is essentially a quark plasma of unbroken symmetry. A second minimum (0 = <j) m in{T)) , 
corresponding to a low temperature hadron phase, occurs at temperatures between Tq and T c , where 



Below To the potential has a local maximum at <j) — and a single minimum at (f> — <p m in{T). 

The resulting differential equations for momentum and energy, neglecting baryon conservation, can be written in 
flux conservative form as 

dtSj + diiSjV*) = -dj\pr - V(4>, T)] - [ K Wd t <p + KWtfd^ + d^V(4>, T)] drf , (2.5) 
d t E + d l {Ev l ) = -\p r -V((f>, T)][d t W + diiWv*)} 

+K{Wd t (j) + Wv t d t 4>f + d^V^, T){Wd t <t> + Wv'd l 4>) , (2.6) 
which, together with the scalar field wave equation, 

d 2 ^ - d l d l( t> = -d^V((f>, T) - KW(d t + v l d t ^) , (2.7) 

and a prescription for the equation of state described below, represent the complete system of equations we solve in 
this paper. In deriving the above equations, we have implicitly assumed a flat space metric, and defined variables so 
that W — u° — — v l Vi is the relativistic boost factor, v l = u l /u° is the transport velocity, S% — Wphui is the 

covariant momentum density, E is the generalized internal energy density 

^ = 3a r T 4 + V{cj>, T) - Td T V(^ T) , (2.8) 
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a r = g*(7r 2 /90) is a radiation constant defining the particle content and degrees of freedom and k is a constant 
coefficient regulating entropy production and kinematic and scalar field energy dissipation. The terms associated with 
k act as an effective friction force at the phase transition surface boundary El . 

The hydrodynamic equations are solved using the Cosmos code [TTL Il2| with time-explicit operator split methods, 
second order spatial finite differencing, and artificial viscosity to capture shocks |l3[ . The scalar field wave equa- 
tion is solved using a second order (in space and time) predictor-corrector integration, sub-cycling within a single 
hydrodynamic cycle since equation ((2.7(1 generally has a tighter restriction on the time step for stability. 

A second formulation of the dynamical equations considered here is based on a simpler conservative hyperbolic 
description of the hydrodynamics equations. In this case, the equations are derived directly by expanding the conser- 
vation of stress-energy T^ v v = into time and space explicit parts. In flat space, this yields the same equations 1(2.5(1 
and 1(2. 7|) for momentum conservation and scalar field evolution, but with 

d t E + d i {£v i ) = -d i [v i {p r -V{(j>, TDj+^ + fcWS^ + W, T)]d t 4> (2.9) 

in place of (|2.6() for the energy variable, defined as 

£ = Wphu - [ Pr - V(4>, T)]. (2.10) 

The Cosmos code uses a non-oscillatory central difference numerical scheme ^4| with second order spatial differencing 
and predictor-corrector time integration with dimensional splitting to solve the hydrodynamics equations in this 
formalism. Artificial viscosity is not needed in this approach. Instead, shocks are captured with MUSCL-type piece- 
wise linear interpolants and high order limiters for flux reconstruction, but without the complexity of Ricmann solvers. 

A final ingredient needed is an equation of state relating p r to energy E and temperature T . Here we assume 
p r = a r T 4 , treating the plasma as a gas of relativistic particles with enthalpy density 

ph = JUL- - Td T V(& T) = Aa r T 4 - Td T V{^ T) = ^-+p r - V{4>, T) , (2.11) 
7 — 1 W 

and adiabatic index 7 = 4/3. The equation of state is then written in terms of computed quantities 

Peff =p r - V(4>, T) = ( 7 - 1) (JL + TdrVtt, T) - V{ff>, T)) - V{ff>, T) , (2.12) 

or 

Peff = Pr - V(4>, T) = (7 - 1) ^ 2 _ (7 _ ^ ) ~ V{<f>, T) , (2.13) 

depending on which energy variable is evolved. The quantity p e // is introduced as an effective pressure for numerical 
convenience in solving equations 1)2. 5[l . 1)2.6(1 and/or 1(2.9(1 . The fluid temperature is computed from either 

E - 3a r T 4 - V((f>, T) + TdrV^, T) = , (2.14) 

or 

£ - V{cf>, T) + TW 2 d T V((t>, T) - 3a r T 4 [jW 2 - (7 - 1)] = . (2.15) 

Equations ((2.14(1 and 1(2.15(1 are solved using a Newton-Raphson method to iterate and converge on the temperature, 
assuming an initial guess of T = (E/SWdr) 1 / 4 . 



III. PERTURBATION ANALYSIS 



In this section we carry out a first-order perturbation expansion of the hydrodynamics equations, and review briefly 
the expected dynamical behavior and stability of phase transitions in the context of linear theory. The general aproach 
follows that presented in references 0,H,B>@- I n addition to elucidating any unstable behavior, the results of this 
section are also used to set up inhomogeneous perturbations as initial data that is evolved numerically in 3IVI 

It is convenient to start with the hydrodynamics equations in conservative hyperbolic form ((2.5(1 and 1(2.9(1 . which 
we rewrite here as 

d t £ + d t {£v l +v l p) = S° , (3.1) 
d t S j +d i (S j v i +S ij p) = E j , (3.2) 
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assuming a constant scalar field potential on either side of the phase transition front, and S° and S- 7 represent any 
source or dissipative terms. Equations l|3.1|) and l|3.2[) are expanded out to first perturbative order assuming two 
dimensional perturbations off a homogeneous background flow of the form 

P = Po + Sp{t,x,y) , (3.3) 
v = v Q i + [Sv x (t,x, y) i + Sv y (t, x.y) j] , (3.4) 

for the fluid pressure and velocity. The corresponding boost factor and conserved variables are, to first order 

VI - v 2 y/l-vl (1 - Vq) 3/2 

t2 



p(-fW 2 -7 + 1) 7p 

1 = 1 = ~P° + 7 TTFI 2T 

7-1 7-1 )(l-«o) 



ZlPoVo 



( 7 -l)(l- Vo 2 ) 2 ' (7-l)(l-« 2 ) 

£> = 9 = t -rr, tT Po y o + Po ov x + v Op + j 0v x , (3.7) 

l-v 2 (7-l)(l-^) V l-«o 



(Jwa + t 7771 2T <^> _ $p , (3.6) 



7P0 x , q m 

5 = i 2 = 7 TT71 5T H > ( 3 - 8 ) 

l-v 2 (7-l)(l-^) 

where ph = (e+p) = "fp/{"f — 1), and «o is the unperturbed fluid velocity. 

Substituting l|3.5|l - (|3.8|l into l|3.1|l and i|3.2[l yields for momentum conservation along the y-axis 



d v^P) + (7 _ 1 ]^_ u g ) + «oftr(^»)] = ^ v , (3-9) 



for momentum conservation along the s-axis 

(7 - 1)(1 - «g)£* = 1VQ d t (6p) + 7Po (l + v 2 )(l - v 2 )- 1 d t (Sv x ) 

+ (7-l+«o) dx(6p) + 2 1Po vo{l - vl)~ l d x {Sv x ) + jp v d y (Sv y ) , (3.10) 

and for energy conservation 

( 7 - 1)(1 - «o 2 )S° = (1 - v 2 + jv 2 ) d t {8p) + 2 1PQ v Q (l - v 2 )- 1 d t (Sv x ) 

+ jvq d x (6p) +7p (l + Vg)(l - Wo)" 1 d x (8v x )+jp d y (Sv y ) , (3-11) 

all to first perturbative order. 

Equations (|3.10(1 and . 1 1|> can be simplified further by first eliminating d y (Sv y ) from the x-momentum equation 

= ( 7 -l)(l-v 2 )(X x -v X°) 

= v dt(Sp) + ph W 2 d t (Sv x ) + d x {S P ) + W 2 ph v d x (Sv x ) , (3.12) 



then eliminating d t (8v x ) from the energy equation 



1 + vi 



= cj(7-l)(l-t$ 

= (l-v 2 c 2 )d t (5p) + vo{l-ci)d x {6p) + ph [d x (Sv x ) + d y (6v y )} . (3.13) 



In writing (|3.12(l and Ij3.13|l we introduced the notation c 2 = 7 — 1 , Wq = 1/(1 — «q) P^o = IPo /(j — 1) , and explicitly 
ignored the source terms E" by setting them to zero. 

The final equations (|3.9|l . 13.12|l . and (|3.13|l can be written conveniently in compact form as 

At&W + A x d x W + A y d y W = , (3.14) 
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where W = (dp, Sv x , Sv y ), and 



A t = 



A x = 



A v = 





ph W 2 
P h W 2 

el) ph c 2 
W$ph v 
P h W 2 v 

ph c 2 s 



1 



(3.15) 



(3.16) 



(3.17) 



Next, following the general procedure outlined in reference we assume a solution of the form ~W(t, x,y) 



w(x)e 



—i(u)t-\-ky) 



with w(x) = ' 



^Rj and the following characteristic equation or dispersion relation for l|3.14|l : 



(w + X*v )[k 2 (l - «g) - c; 2 (w + X*v f + (A* + wu ) 2 ] = . 

Equation l|3.18|l has either three distinct roots which agree with those found by 0, || 

u 
v 

1 



or two roots 



(c s - l)u w ± c s (l - v )Jlu 2 + 



a; = — , 



w fc 2 w(l + u 2 ) 



2a; 



2v 




if «o = c a . The corresponding eigenvectors of 13.14fl are 



Ri 



R 



R. 



/ (ph a ujc 2 + ph v c s n) / (k(v 2 - c 2 s )) 

I ( WoW + c s n) (i - Wo 2 )/(fc(c 2 - v 2 )) 

(photic 2 - phoVoc s ft) /(k(v 2 - c 2 )) 
(v u - c s n) (1 - W 2 )/(fc(c 2 - v 2 )) 



1 



(3.18) 

(3.19) 
(3.20) 

(3.21) 
(3.22) 

(3.23) 
(3.24) 
(3.25) 



where we define fl = [lo 2 + k 2 (v 2 - c 2 )/(l - v^)] 1 / 2 in equations l|3.24|l and Q3.25fl . 

To determine whether any unstable modes of W exist and are consistent with imposed boundedness conditions 
at x — > ±oo, it is convenient to construct a coordinate system centered and moving with the frame of the interface 
between quark and hadron phases at x = 0. Assuming the high temperature quark phase is to the left of the interface 
(x < 0), and the hadron phase to the right (x > 0), we can easily determine if any unstable modes (defined by 
Im(w) > 0) obey the conditions W — > as x — > ±oo. As concluded in references 0,1^], if Im(w) > in the quark 
region of a detonation front, then Im(A*) < and therefore a,j = for all j in order for the solutions to be bounded at 
x — > — oo. Hence unstable modes do not exist in regions ahead of nucleating bubbles separated from the quark phase 
by a detonation front. However, in the hadron region to the right (x > 0) of a detonation front, the boundedness 
conditions at x — > +oo do not rule out unstable modes to linear order. Also, in the case of deflagration fronts, unstable 
modes are realizable on both the quark and hadron sides, at least for cases in which surface tension, dissipation, and 
strong thermal coupling effects are negligible Q, H, Q . It is not entirely clear what role these effects play in stabilizing 
the transition, and certainly higher order nonlinear effects are not accounted for in this analysis. Table [I] summarizes 
which of the regions and hydrodynamic states may potentially give rise to unstable modes with Im(w) > 0, at least 
to linear order. 
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TABLE I: Summary of parameters yielding bound solutions that are susceptible to unstable modes. The entries represent a 
general set of conditions that allow unstable modes (Im(w) > 0) to first perturbative order for each of the phases and front 

types. These conditions are used to construct initial data for the multi-dimensional numerical stability studies presented in 
section tllVI 

Quark region (x < 0) Hadron region (x > 0) 

detonations: detonations: 

Vq > C B V h < Vq 

ai = A* = —w/vo 

a- = a- = 

a+ = MA+) < 

deflagrations: deflagrations: 

Vq < C s V h > Vq 

ai = A* = —cj/vq 

Im(A*)>0 a_=0 

a+ = Im(A+) < 



IV. NUMERICAL RESULTS 
A. Initial Data 

The initial data has just one free dimension, energy, which we use to specify the critical transition temperature 
T c and normalize it to unity. It is also convenient to define a characteristic length scale as the radius of a spherical 
nucleated hadron bubble in approximate equilibrium with an exterior quark plasma. This is estimated by balancing 
pressure forces acting from the quark and hadron sides and including the effect of surface tension, resulting in 
R eq — 2a/(ph — p q ) for the equilibrium radius, where a is the surface tension, and ph and p q are the pressure 
fields from the hadron and quark sides respectively. For the preliminary one-dimensional calculations in iflVBI the 
computational box sizes are set to a multiple of the equilibrium bubble radius L = xR eq , with x generally in the range 
20 - 2000. Typical length scales or box sizes in these calculations vary from a few hundred to a few thousand fermi. 
By comparison, the proper horizon size at the time of the QCD phase transition is 

r* i6 

d H = 2t age = R(t) / dt'/R{t') = - - km , (4.1) 

Jo v^/51.25 (TM e wl50 MeV)^ 

or 16 km for = 51.25 and TmcV — 150 MeV. In writing (|4.1|l we have assumed R(t) is the cosmological scale factor 
in the isotropic FLRW background model dominated by radiation energy density of the form e = a r T 4 = eoR~ 2 such 
that 



tage ~ V 32irGa r ~ ^g,/5l.25 (T Me y/150MeV) 2 " 

is the age of the Universe as a function of temperature. 

We consider two background temperatures for the super-cooled initial state triggering the shock and phase front 
propagation: T = 0.9T C and T = 0.9943T C . The latter temperature is chosen to match the 1+1 dimensional cases 
considered in reference 0, and provides a useful benchmark against which some of our results can be compared. 

To keep the number of numerical simulations down to a reasonable level, we fix the following additional parameters 
for all calculations and for both temperature cases: critical transition temperature T c = 150 MeV corresponding to 
the QCD symmetry breaking energy occurring when the Universe was approximately t age = 2.7 x 10 -5 s old (see 
equation l|4.2|l ): particle degrees of freedom = 51.25 to be consistent with previous studies [l^; surface tension 
fj = 0.1 T c ; correlation length l c = 1.0/T C ; and latent heat L = 2 Note that for both temperatures considered 
here, the interaction potential (equation (|2.2|l and Figure ^1 has two local minima (at <j> = and <f> = (j) m i n (T) > 0, 
which allows for the quark (cold phase with <f> = 0) and hadron (hot phase with cf> = <fi m in(T)) states to co-exist. 

The remaining parameter, the dissipation constant k, is varied over the range 0.1 < n < 10 in the one-dimensional 
calculations to determine the regimes in which detonations and deflagrations are triggered, and to compute the 
stability parameters from linear theory. 
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The scalar field <j> is initialized according to 



= Ml + A)— I j^-^ , (4.3) 



with 



= i-AcothV/?(r2-r,?) x/(i + A)) 

l-coth 2 (V/3(T 2 -T 2 ) + A)) ' 

A = 9A/?(T 2 — Tq) / '(2a 2 T 2 ) , and x = L is the right-most edge along the x-axis representing the computational box 
size. These expressions define a hadron phase with </> > at L — x — > 0, and quark phase with — > as L — a; ^ oo, 
with an exponentially steep interface boundary of width set by the model parameters. The numerical constants <p c 
and A = 0.05 are used to perturb the solution out of dynamical equilibrium by varying both the amplitude and width 
of the field in the hadron phase. For most of the calculations we set <p c — 1. However, for the larger k cases we found 
it useful to specify <p c such that <p(x = L) w (j) min (T), in order to reduce the transient interval between the initial 
configuration and the time that the phase front fully develops. 



B. One-Dimensional Results 



The primary motivations for the one-dimensional calculations presented here were to narrow the parameter space 
and help illuminate the two-dimensional studies presented in the following section for the different background tem- 
peratures considered. In cases where the temperature was initialized to T = 0.9T C , we used anywhere from 2 x 10 3 
to 3.2 x 10 4 zones to resolve a length scale of L = 2000i? e(? = 2190 fm. For runs in which T — 0.9943T C we used from 
5 x 10 2 to 4 x 10 3 zones to resolve L — 20R eq = 452 fm. In these ID parameter studies, we also considered an initial 
background temperature of T = 0.95T C . These runs were useful for characterizing some of the intermediate states in 
the 2D runs below. For these T — 0.95T C runs we used from 1 x 10 3 to 4 x 10 3 zones to resolve L = 250R eq = 1157 fm. 
The higher resolutions in each temperature case were generally required for the larger values of n in order to maintain 
reasonable zone coverage in the hadron phase, since the bubble growth velocity was sometimes substantially less than 
the shock front velocity, or near the transition from a detonation front to a deflagration front, as detailed below. For 
this first group of ID runs, the code was allowed to run for approximately one sound-crossing time t s = L/c s , where 

c s = 1/V3. 

Figure |21 shows the phase front or bubble growth velocity normalized to the sound speed Vf/c s , and the quantity 
r\ = —T c Vfdvf/dT q as a function of the dissipation parameter k for temperatures T/T c = 0.9, 0.95, and 0.9943. 
According to reference Q, r\ determines the stability of bubbles and we compute it using the approximation r\ = v'j/v 2 , 
where w 2 = (1/2)(1 — T 2 /T 2 ). Notice the T = Q.9T C curve switches from a detonation to a deflagration at k f» 1 
where Vf/c s crosses unity, and the T = 0.95T C curve switches at n w 0.3. However, the precise transition points 
are highly sensitive to the resolution and duration of the simulations, and what might appear as a deflagration at 
one resolution may turn out to be a detonation at another. This is attributed to the solutions approaching either a 
Jouget or temporary strong detonation state (which eventually decays into a weak deflagration) as the parameter k 
is increased. In particular, the post-shock velocity for the T = 0.9T C case computed over short time intervals in the 
frame of the front decreases from being supersonic at low values of k, to about sonic speed at K = 0.875, then subsonic 
at higher values but less than unity (e.g., k — 0.9), corresponding to weak, nearly Jouget, and strong detonations, 
respectively. Solutions for which n is close to unity thus represent a regime in parameter space that does not allow 
stable detonation modes since strong detonations are forbidden modes of propagation (see, for example, 01)- We 
note that the conservative hyperbolic form of the hydrodynamics equations and the non-oscillatory central difference 
methods for solving the equations generally out-perform artificial viscosity methods in resolving and maintaining 
stable evolutions near the critical Jouget state. 

Table ITT1 summarizes our results for the potentially unstable, rj < 1, deflagration runs. Th and T q in Table UTI are 
the average temperatures of the hadron and quark phases, respectively; and Vh and v q are the hadron and quark 
velocities on either side of the phase front as measured in the rest frame of the moving front. The critical wave 
number, k c = (fi — l)ph q v q /a, defines the limit above which the system is stabilized by surface tension where 
A* = Vh/v q and ph q is the enthalpy density of the quark phase. Modes with wavelengths A < A c = fc" 1 are thus 
expected to be stable. The instability growth timescale, r, is estimated as Q 

k Q v q , (4.5) 



1 

A* 

A* 



1 + fi ak 
H ph q v 2 q 



1/2 
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FIG. 2: Bubble growth velocity normalized to the sound speed Vf/c s , and stability parameter tjasa function of the dissipation 
constant k for temperatures T/T c = 0.9, 0.95, and 0.9943. The T = 0.9T C (T = 0.95T C ) curve switches from a detonation to a 
deflagration at n « 1 (k ~ 0.3). 



TABLE II: Hadron temperature Th, quark temperature T q , front velocity Vf , hadron velocity relative to front Vh, quark velocity 
relative to front v q , critical stability wavenumber fc c , fastest growth wavenumber fco, linear perturbation growth time r, and 
stability parameter r\ as a function of dissipation constant k for initial temperatures T/T c = 0.9, 0.95, and 0.9943. The results 
are shown only for the potentially unstable ID deflagration runs in Figure [5] that satisfy r\ < 1. 



T/T c 




T h /T c 


T q /T c 


\Vf/c s \ 


Vh 


v q 


h c 


fc 


T 




0.9 


5.0 


0.90639 


0.92971 


0.431 


0.249 


0.193 


1.816 


0.933 


8.25 x 10 1 


0.911 




6.0 


0.90378 


0.92449 


0.364 


0.210 


0.164 


1.250 


0.642 


1.46 x 10 2 


0.608 




7.0 


0.90275 


0.92264 


0.339 


0.196 


0.153 


1.071 


0.550 


1.83 x 10 2 


0.516 




10.0 


0.89914 


0.91717 


0.260 


0.150 


0.118 


0.613 


0.314 


4.17 x 10 2 


0.285 


0.95 


5.0 


0.95436 


0.96233 


0.258 


0.149 


0.127 


0.546 


0.278 


6.81 x 10 2 


0.602 




7.0 


0.95246 


0.95964 


0.203 


0.117 


0.100 


0.331 


0.170 


1.43 x 10 3 


0.348 




10.0 


0.95072 


0.95747 


0.158 


0.091 


0.077 


0.198 


0.101 


3.09 x 10 3 


0.199 


0.9943 


1.5 


0.99633 


0.99686 


0.082 


0.048 


0.043 


0.043 


0.022 


4.33 x 10 4 


0.722 




2.0 


0.99596 


0.99646 


0.071 


0.041 


0.037 


0.031 


0.016 


7.06 x 10 4 


0.479 




3.5 


0.99527 


0.99576 


0.043 


0.025 


0.022 


0.012 


0.006 


2.66 x 10 5 


0.145 




5.0 


0.99491 


0.99541 


0.033 


0.019 


0.016 


0.007 


0.004 


6.82 x 10 5 


0.077 




10.0 


0.99439 


0.99492 


0.016 


0.009 


0.007 


0.002 


0.001 


5.95 x 10 6 


0.017 



where 



ka = - 



9 V 1 + It 



3/x 



3/i — — 



1/2' 



PK V q 



(4.6) 



is the wave number of the mode with the shortest growth time scale. 

Table lllll summarizes the results computed for the detonation cases with temperature T = 0.9T C . Here Vf represents 
both the front and shock velocities, and phh is the enthalpy density of the hadrons immediately behind the phase 
front. The quantity cr/ph has dimensions of length, and provides a natural scale that gauges the relative importance 
of surface tension. It is also applied as a dimensional scaling variable in the calculations of reference 0, which we 
use as a guide to estimate perturbation growth rates and physical run times for the two-dimensional calculations 
presented in fllVCI 
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TABLE III: Front velocity Vf and enthalpy to surface tension ratio (phh)/<J on the hadron side as a function of dissipation 
constant k for a selected set of ID detonations with initial temperature T = 0.97^ 



T/T c 
0.9 



0.1 

0.3 

0.5 

0.75 

0.8 



1.85 
1.71 
1.57 
1.36 
1.33 



{phh)/cr 

191.7 

195.7 

201.1 

216.2 

220.4 



In Figures 13151 we present more detailed results of three particular ID runs that will help illuminate the 2D 
calculations with similar parameters presented in the next section. Figure [3] shows a series of outputs of the scalar 
field and temperature for a weak (subsonic post-front velocity relative to the phase front) deflagration case with 
T = 0.9T C and n = 7. This run used 2000 zones to resolve a length scale of L = 1095 fm, and was allowed to run 
for slightly longer than two sound-crossing times so that the shock leading the deflagration front propagates across 
the grid twice. These tiled displays show a leading shock front (in the temperature graphs) originating from the 
right and traveling to the left with velocity slightly higher than sound speed v s = 1.03c s , heating up the quarks just 
ahead of the cool hadron phase (first row) . The deflagration front (separating the two regions with different scalar 
field values) travels to the left at a speed slower than the leading shock front, Vf/c s ~ 0.34. Due to the reflection 
boundary conditions imposed on all our calculations, a mirrored shock collision occurs when the leading shock hits 
the left boundary (second row). The reflected shock then travels to the right, eventually colliding with the phase front 
moving in the opposite direction (third row). The shock heats up the fluid and reduces the scalar field as it passes 
through the hadron phase. Upon impact with the phase front, the shock/front interaction also generates a rarefaction 
wave traveling to the left which cools off the quarks in its path. Finally the hadron matter at the right end is heated 
further as the reflected shock hits the right boundary and collides with another incoming shock (fourth row). 

Figure 0] shows a similar series of displays for a weak (supersonic post-shock velocity relative to the phase front) 
detonation case with T = 0.9T C and K = 0.5, run for almost two sound-crossing times. This run used 2000 zones to 
resolve a length scale of L = 1095 fm. The outputs begin in the first row with an early undisturbed state showing the 
leading shock/detonation front originating at the right end of the grid and traveling to the left at speed Vf/c s w 1.6. 
A rarefaction wave follows the shock and separates two hadron regions at two slightly different temperatures. The 
second row shows the shock-shock collision which occurs when the leading shock front reaches and reflects off the left 
boundary. This interaction heats up the fluid to temperatures higher than T c , and generates a quark region at the left 
boundary. As the reflected shock travels to the right, the quark region grows. The third row shows the reflected shock 
passing through the rarefaction wave separating the two hadron states. This rarefaction wave cools the newly formed 
quark region, converting quarks back into hadrons. The fourth row shows the interaction of reflected rarefaction 
waves at the left boundary. The complete reflected state showing the thermal distortion and shock profiles after the 
entire domain is converted to the hadron phase is displayed in the fifth row. This compact structure is allowed to 
traverse the grid and interact through each sequence again at the right boundary (sixth row). The complete reflected 
state from the right boundary is shown in the final row. 

Figure [S] shows outputs from two phase front systems originating in separate regions at different temperatures and 
are thus out of initial thermal equilibrium: a deflagration front at T — 0.9943T C at the left end, and a detonation 
front at T = 0.9T C at the right end. This run used 2000 zones to resolve a length scale of L = 1806.1 fm. The 
dissipation constant is k — 0.5 across both regions, and the temperature discontinuity is initiated at the center of the 
grid. The right end is initialized to the same state as the previous detonation case in Figure 0] The left end is similar 
to the weak deflagration case of Figure |3 but with a smaller dissipation constant and a higher initial temperature. 
The first row illustrates an early stage where the deflagration front has formed at the left end with Vf/c s w 0.14, 
the shock/detonation front and rarefaction wave have formed at the right end with v s = Vf ss 1.6c s , and a thermal 
shock front and rarefaction wave are traveling in opposite directions from the central thermal discontinuity at about 
the sound speed. The second row shows the collision of the detonation and thermal shocks. The third row shows the 
passage of the thermal shock through the detonation rarefaction, and the collision of the deflagration shock with the 
thermal rarefaction wave. The fourth row shows the interaction of the deflagration front with the thermal rarefaction 
at the left end, and the collision of the deflagration and detonation shocks at the middle. As the rarefaction wave 
passes through the hadron phase behind the deflagration front at the left end, the phase transition is accelerated to 
near sound speed. The fifth row shows the accelerated shock and deflagration front features traveling to the right at 
velocities v s 0.67 = 1.16c s and vj « 0.51 = 0.88c s . Notice that the accelerated front velocity is close to, but slightly 
faster than the deflagration velocity predicted by the T — 0.95T C curve of Figure|3 which yields Vf — 0.74c s = 0.43. 
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FIG. 3: Scalar field (left column) and temperature (right column) as a function of position for a ID deflagration case with 
T — 0.9T C and k = 7. Top row: Early undisturbed deflagration and shock fronts at t = 3.2 x 1Q -21 s. Second row: Shock 
collision in the quark phase at the left boundary at t = 6.6 x 10 -21 s. Third row: Reflected shock interacting with deflagration 
front at t = 1.1 x 10 -20 s. Fourth row: Interaction of two reflected shocks in the hadron phase at the right boundary at 
t = 1.4 X 10" 20 s. 



The sixth row shows the interaction of the newly accelerated front with the detonation propagating from the right 
end. This interaction heats up and decomposes the hadrons back into quarks over a substantially large region at the 
center. The newly formed quark region gradually cools with the shock passage and rarefaction wave interactions. 
Profiles and distributions of quark and hadron regions become increasingly complex as indicated in the final row due 
to the various shock, phase, and sound wave interactions as they reflect off the boundaries and propagate through the 
grid. 

We point out that, in general, Figure |2 cannot always be used to predict reliably the accelerated velocity or mode 
of propagation from rarefaction wave/front interactions. For example, consider the same initial configuration used in 
producing Figure but with friction parameter k = 0.2. In this case Figure [3 predicts, for a mean field temperature 
T = 0.95T C , a supersonic front velocity off/ — 1.56c s . However, the actual (numerical) accelerated solution remains a 
weak deflagration with subsonic velocity Vf = 0.54 = 0.94c s , only slightly faster than for the k = 0.5 case. This holds 
even in the limit k — > 0: the accelerated front velocity approaches but does not exceed the sound speed as k — > 0, and 
the propagation mode remains a deflagration with clearly separated phase and shock front features. 
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FIG. 4: Scalar field (left column) and temperature (right column) as a function of position for a ID detonation case with 
T — 0.9T C and k = 0.5. Top row: Early undisturbed detonation front plus rarefaction at t = 2.2 x 1CF 21 s. Second row: 
Shock-shock interaction at left boundary at t — 4.8 x 10~ 21 s. Third row: Reflected shock passing through the rarefaction wave 
at t = 5.7 x 10~ 21 s. Fourth row: Interaction of the rarefaction waves at left boundary at t = 6.6 x 10~ 21 s. Fifth row: Reflected 
state from left boundary at t = 9.7 x 1CF 21 s. Sixth row: Interaction of reflected states at right boundary at t = 1.2 x 1CT 20 s. 
Seventh row: Reflected state from right boundary at t = 1.5 x 1CT 20 s. 
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FIG. 5: Scalar field (left column) and temperature (right column) as a function of position for a ID interaction of shocks 
and phase fronts originating in two separate regions at different temperatures: a deflagration system at T = 0.9943T C at the 
left end, and a detonation configuration at T — 0.9T C at the right end. The dissipation constant is k = 0.5 everywhere, and 
the temperature discontinuity is initiated at the grid center at t = 0. Top row: Early undisturbed deflagration and shock 
fronts (left), detonation front plus rarefaction (right), and thermal shock front plus rarefaction (center) at t = 1.8 x 10 -21 s. 
Second row: Interaction of detonation front and thermal shock at t = 2.2 x 10 -21 s. Third row: Interaction of the detonation 
rarefaction and thermal shock, and thermal rarefaction and deflagration shock at t = 3.1 x 10 -21 s. Fourth row: Interaction 
of deflagration front and thermal rarefaction, and deflagration shock and detonation shock at t = 4.8 x 10 -21 s. Fifth row: 
Acceleration of deflagration front by thermal rarefaction at t = 5.3 x 1CP 1 s. Sixth row: Interaction of newly accelerated front 
and detonation front at t = 7.5 x 1CF 21 s. Seventh row: Interactions of multiple shocks resulting in two separate quark regions 
at t = 1.1 x 10 -20 s. 
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TABLE IV: Summary of parameters for the two-dimensional runs. 



Run 


T IT 

J- / J-c 




\ 

A c 




I -.viri Tlimoncinnc 
Vjrl 1L1 .Lx 11 lie 11 OIL* Ufa 




^Itnn Timp 








(fm) 


00 


(fm) 


(zones) 


00 


A 


0.9 


7.0 


1.2 


1.3 x 10" 21 


278.03 x 6.14 


1280 x 64 


3.2 x 10~ 21 


B 


0.9943 


1.5 


31. 


2.6 x 10" 19 


3059.83 x 152.99 


1280 x 64 


7.1 x 10 _2 ° 


C 


0.9 


0.5 


6.6 


2.2 x 10~ 22 


302.34 x 32.71 


1024 x 256 


3.5 x 10~ 21 


D 


0.9 


7.0 


1.2 


1.3 x 10~ 21 


1000 x 1000 


512 x 512 


1.2 x 10~ 20 


E 


0.9 


0.5 


6.6 


2.2 x 10 -22 


1000 x 1000 


1024 x 1024 


1.2 x 10~ 20 


F 


0.9943/0.9 


0.5 






1806.1 x 451.53 


2048 x 512 


2.1 x 10~ 20 



C. Two-Dimensional Results 

We focus here on extending to two dimensions several of the more interesting parameter combinations presented 
in section i ilVBI that may potentially give rise to unstable behavior as predicted by linear theory. In particular, 
we consider six calculations: runs A and B are single deflagrating fronts; run C is a single detonation transition; 
runs D and E are interactions between a planar and smaller spherical nucleating bubble of deflagration (run D) and 
detonation (run E) types; and run F is the interaction of a deflagration system with a detonation front nucleating 
from two regions out of thermal equilibrium. Table ITVl summarizes the parameters used in each of these runs. 

For all of these simulations we initialize the data with various perturbations included. First, the planar fronts in 
each problem are perturbed with a sinusoid of wavelength A rs 5A C , where A c is the critical wavelength calculated from 
our ID studies. The amplitude of this perturbation is Sx/x — 0.05, and the wavelength A is used to set the physical 
size of the grid parallel to the phase front. 

Transverse and longitudinal inhomogeneous fluctuations are also imposed on all the initial data using the per- 
turbation solutions discussed in flllll In particular, E — Eq + SE, v x = Sv x , and v y = Sv y , with Eq — 
3a r T 4 + V(cj), T) - Td T V(cj), T), where 

/ (<y-l)8(E/W)\ ^ 
Sv x \ =e- i ^ t+ky ^J2 a ^i e (X1)Xe ( ^ )X ' ( 4J ) 

V SVy J * 

A* are the eigenvalues, and Rj = Rj/|Rj| are the normalized eigenvectors. The expansion coefficients are defined as 
dj — min(Aw , AEq) with amplitude constant A = 0.1 (if they are not set to zero because of boundedness constraints 
as discussed in Mlllf) . vq is the background average velocity in either the quark or hadron regions as defined in £ 11111 
but implied here to be the velocity component orthogonal to the interface boundary and measured in the rest frame 
of the unperturbed surface. We take vq = c s + 0.5(1 — c s ) or vq — 0.1c s for the unperturbed velocity in the detonation 
or deflagration cases, respectively. Finally, we also impose random fluctuations in the background temperature with 
peak excursions of ST/T = 0.005. 

1. Single Front Perturbations 

The growth and/or decay of perturbations are tracked by calculating the power spectral density (PSD) of various 
wavelength modes as a function of time. Here the PSD is defined as P(k) = 2\Ck\ 2 /N 2 , where N is the number of 
zones parallel to the front, k = 1, 2, • • • , (N/2 — 1) is the wavenumber, and Ck is the discrete Fourier transform. We 
calculate the PSD for the phase front position as well as the integrated scalar field and integrated temperature behind 
the front as a function of transverse coordinate (along the y-axis). 

Considering the results of section fllVBI we choose for the first deflagration case (run A) the parameters T = 0.9T C 
and k = 7, which results in a growth timescale r that is computationally realizable. We use this growth time to set 
the run time (t w 2.5r w 3.2 x 10 -21 s) and the physical length scale (~ c s r) perpendicular to the phase front. This 
problem is resolved with 64 zones parallel to the front and 1280 zones perpendicular to the front, providing similar 
grid spacing in each of the two dimensions. Figure shows the power spectral density as a function of time for the 
k = 1 and k = 2 modes of the front position and the integrated scalar field behind the front. The data for the phase 
front is terminated at t w 0.03r «5x 10~ 23 s (or t — 12 in code units), at which time the perturbations in the front 
position are too small to be resolved at the grid resolution scale. 
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FIG. 6: Power spectral density as a function of time for the k — 1 and k = 2 modes of the front position and integrated scalar 
field behind the front for the T = 0.9T C , k — 7 deflagration case (run A). 




FIG. 7: Same as Figure^] except for the T = 0.9943T C , « = 1.5 deflagration case (run B). 



We also consider a second deflagration case with T — 0.9943T C and k = 1.5 (run B) for comparison, despite the fact 
that the predicted growth timescale is well beyond what we can simulate in real time. Instead we set the physical 
length scale perpendicular to the phase front to be twenty times the scale parallel to the front, which is « 5A C as 
before. We then fixed the run time to be t w 4t s w 0.25r w 7.1 x 10 -20 s, where t s — L/c s w 100A c /c s is the 
sound-crossing time perpendicular to the front. This physical time scale gives us a reasonable run time to formulate 
conclusions about the general behavior of perturbations across the shock and phase fronts. Figure shows the PSD 
as a function of time for this simulation. Data for the phase front is terminated at t as 0.04r as 1CP 20 s (t = 2518 in 
code units) once the perturbations in the front position are smaller than the grid resolution. Note that the physical 
scale perpendicular to the front for this run is about a factor of ten larger than for runs A and C. This implies that 
the initial perturbation in the phase front is also a factor of ten larger, accounting for the greater power in each of the 
modes and partially explaining the longer decay time. The decay time is also influenced by the slower front velocity 
and longer wavelength in this case. 

For the detonation case (run C), we choose T = 0.9T C , k = 0.5 and a wavenumber k = 0.001 (ph^/a), which, 
according to Figure 2 in gives a growth time of r ~ 1000 (a/phh) for Chapman- Jouget detonations. Using the 
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tabulated value of phh/cr in Table ITTT1 we set the simulation run time to 16r = 3.5 x 10~ 21 s. Figure [S] shows the PSD 
as a function of time. Again, data for the planar front position is terminated atf^T (t — 45.4 in code units) when 
the perturbations are too small to be resolved within a single grid cell. 

Although the shape function of the phase front and the individual modes in each simulation are observed to oscillate 
in time, the oscillations are strongly damped over a period that is significantly shorter than the characteristic growth 
time predicted by linear theory. We thus find no evidence of unstable behavior in any of the perturbed planar 
deflagration or detonation cases we have investigated, a result that is generally consistent with the linear analysis and 
conclusions of Huet et al Ijl . 



2. Multiple Front Interactions 

Here we are interested in the question of whether mixing can occur through a turbulence-type mechanism triggered 
by shock or pressure wave distortions across phase boundaries. For this, we first consider two-dimensional interactions 
of multiple nucleating regions of varying sizes and radius of curvature. Runs D and E look at the interaction of a 
perturbed planar front system with a smaller nucleating bubble region. Run D is a deflagration system with T — 0.9T C , 
k = 7, and run E is a detonation system with T = 0.9T" C , k = 0.5. The physical dimensions of the grid used in both 
problems is 1000 x 1000 fm, and is resolved with 512 x 512 cells in run D and 1024 x 1024 cells in run E. In both runs, 
the small nucleating bubble region has an initial radius of 30 fm, and the simulations were run to t w 2t s , where t s is 
the sound crossing time across the grid. 

Figures l9l and ITUI show the scalar field and the temperature for run D, up to time t ss 1.8t s in the final image. 
These figures show behavior similar to the one-dimensional plots of Figure |3J despite the greater complexity of the 
shock and phase geometries, and the additional multi-dimensional perturbations. The separation of phase and shock 
fronts is clearly evident in the temperature contours, as are the reflected shock and rarefaction waves generated in the 
collisions. Perturbations along and behind the phase fronts are observed to decay in time as expected from the results 
of section [TV C II and we observe no signs of turbulent behavior. In this case, the hadron phases merge smoothly as 
the bubbles grow, leaving behind only a spectrum of thermal fluctuations from the initial perturbations, geometrical 
scales, oblique shock interactions, and an increasing mean temperature with each global shock passage. This scenario 
is thus consistent with the standard picture of the QCD transition in which hadron bubbles nucleate and expand as 
spherical subsonic condensation discontinuities undergoing a process of collisions, regular (non-turbulent) coalescence, 
and quark droplet decay, to eventually fill the universe. 

Figures |lll and |12I show the scalar field and temperature for run E up to time t « 1.6t s in the final image. In 
contrast to the deflagration results of Figures 151 and ITUI these images exhibit more complex behavior. In particular, 
front collisions result in the continuous generation of coherent quark "nuggets" formed by the greater entropy heating 
of hadrons at shock contact regions, a result consistent with the one-dimensional calculations shown in Figure The 
nuggets are generally formed over spatial scales set by the length of the thermal plateau separating the leading shock 
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FIG. 9: Contour plots of the scalar field for an interacting perturbed plane and a smaller nucleating bubble in the deflagration 
limit with T = 0.9T C and k = 7.0. Results are shown at times 0, 2.0, 4.1, 6.1, 8.1, and 10.2 x 10 -21 s for grid dimensions of 
1000 fm along both the x and y axes. All images are scaled to the same global color map scheme. 

from the rarefaction wave (or equivalently by the coherence length of the hot leading phase of hadrons at the time of 
impact). This, in turn, is set by the mean free path between hadron bubbles or thermal discontinuities and the strength 
of the detonation front which sets the relative velocity difference between the shock and rarefaction fan. A simple 
dimensional argument suggests that maximum quark fragment sizes are expected to be of order 61 = SvSt ~ L/6, or 
roughly 16% of the box size. In this crude estimate we assume for the velocity difference determining the plateau 
width Sv ~ Cs/3 as found in the ID numerical solutions, the time between collisions 8t ~ Lj (2c s ), and L/2 represents 
the free path length between bubbles in this simulation. The numerical results find quark regions ranging in size 
from the smallest resolved scale of approximately one fermi determined by the cell size, up to a few hundred fermi, 
consistent with the above expected result. Once formed, the nuggets eventually decay away over roughly a sound 
crossing time due to adiabatic cooling and the generation of rarefaction waves at phase boundaries. Nugget shapes are 
generally determined by local shock and phase front geometries, but invariably they tend to become more spherical 
as the droplets decay and surface tension effects become significant enough to help erase surface perturbations. 

A final calculation presented here is run F, which considers the interaction of a deflagrating system with a detonation 
front nucleating from two different regions in space out of thermal equilibrium. The parameters are the same as those 
used in the one-dimensional calculations of Figure[5]in scction lTV Bl but we also introduce a phase shift of ir/2 between 
the two fronts along the transverse direction. Equal grid spacing is used in the two dimensions here, but with four 
times as many zones in the direction perpendicular to the fronts. The physical size of the grid is set to 1806.1 x 451.53 
fm and resolved by 2048 x 512 zones, chosen as a reasonable balance between resolution and run time. The physical 
run time of the final displayed image is t = 1.58t s = 1.65 x 10~ 20 s. The images in Figures H~3l - 1141 for this case show 
many of the same features exhibited by the interaction of two detonation fronts, in particular the decomposition of 
condensed hadrons into quark droplets or decaying nuggets at shock collisions. The spontaneously generated nugget 
regions remain (or become) regular during their lifetime, as do the phase boundaries of the larger initial phase fronts. 

Another interesting feature found in this case is the deflagration 'instability' also noted in the one-dimensional 
results of Figure [5] Although precise quantification of thermal and kinematic states ahead and behind the accelerated 



17 




0.33 0.90 0.92 0.94 0.9 6 0.93 1.00 



FIG. 10: As Figure^] but showing the fluid temperature in units of the critical temperature (T c — 1) at the same corresponding 
times. 

phase front is difficult due to the multi-dimensional perturbations, the phase front is somewhat easier to characterize 
since perturbations are not as pronounced in <f>. The deflagration front at the left side, as measured halfway along 
the y-axis, is accelerated from an initial velocity of vj = 0.35c s = 0.2 to a supersonic speed of about 1.2c s = 0.68 
by the passage of a rarefaction wave generated at the thermal discontinuity. However, the accelerated front velocities 
measured along the top or bottom of the y-eods is smaller and closer to the one-dimensional result described in section 
II V Bl The midpoint of the phase front, which lags behind the mean front position as shown in Figure 1131 thus 
undergoes a greater acceleration due to surface tension and mode dissipation effects. Also, transverse perturbations 
affect the flow of the hadron phase by 'folding' the cold phase into itself, creating hadron domains ahead of the main 
front which enclose decaying quark regions and eventually merge to effectively increase the front propagation velocity. 

V. SUMMARY 

We have performed both one- and two-dimensional numerical simulations of first order quark-hadron phase tran- 
sitions in the early Universe. The primary goal of these studies was to explore the nature of the phase transitions 
beyond linear stability analysis, and determine if the interface regions between phases are inheritly stable or unstable 
when the full nonlinearities of the relativistic scalar field and hydrodynamic system of equations are accounted for. 
Whether interface boundaries are stable or unstable can have important consequences in our understanding of the 
standard picture of phase transitions and for the evolving baryon number density distribution, which we assume 
in this paper evolves hydrodynamically and coupled thermally and kinematically to a nonlinear scalar field with a 
quartic self-interaction potential. We used results from linear perturbation theory to define initial fluctuations on 
either side of the single phase front simulations, and evolved them numerically in time for both weak deflagrations 
and detonations. No evidence of unstable behavior is found in any of the isolated, single perturbed planar front cases 
we considered, despite the fact that the cases we chose to study are predicted to be unstable according to several 
linear analysis calculations. Indeed, the power spectra computed for the phase surface boundary and integrated fields 
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FIG. 11: Contour plots of the scalar field for an interacting perturbed plane and a smaller nucleating bubble in the detonation 
limit with T = 0.9T C and n = 0.5. Results are shown at times 0, 1.4, 2.7, 4.1, 5.4, 6.1, 7.1, 8.1, and 9.1 x 1CT 21 s, on a grid of 
dimension 1000 fm in each direction. 



19 




0.33 0.92 0.95 0.99 1.03 1.06 1.10 

FIG. 12: As Figure fTTI but showing the temperature. 
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FIG. 13: Contour plots of the scalar field for interacting perturbed plane fronts originating in two separate regions at different 
temperatures: a deflagration front at T — 0.9943T C at the left end, and a detonation front propagating from the right end at 
T — 0.9T C . The dissipation constant is k — 0.5 everywhere on the grid of dimensions 1806.1x451.53 fm, and the temperature 
discontinuity is initiated at the grid center along the :r-axis. Results are shown at times 0, 4.0, 8.0, 12.5, and 16.5 x 10~ 21 s. 
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FIG. 14: As Figure ED but for the fluid temperature. 
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behind the phase front all suggest that perturbations decay fairly rapidly, in fact much shorter then the growth time 
predicted by perturbation theory. 

We are also motivated to determine whether phase mixing can occur through a turbulence- type mechanism triggered 
by shock proximity or the disruption of interfaces by pressure or shock waves. To investigate this scenario, we 
considered three cases: the first two involve the interaction of perturbed planar fronts with smaller nucleating hadron 
bubbles, one being entirely a deflagration scenario and the other a detonation; and a third case simulates the interaction 
between deflagration and detonation systems arising from two regions of space which have super-cooled to different 
temperatures and are thus out of initial thermal equilibrium. In the first case of planar-bubble deflagrations, the 
results are consistent with the standard picture of cosmological phase transitions in which hadron bubbles expand 
as spherical condensation fronts and undergo a process of regular (non-turbulent) coalescence, eventually leading 
to collapsing spherical quark droplets in a medium of hadrons. This behavior is generally supported also by the 
second case considering the interaction of multiple detonation bubbles. Although the evolutions in this second case 
are complicated by greater entropy heating from shock interactions, which contributes to the irregular destruction 
of hadrons and the creation of quark nuggets, the interfaces between phases remain coherent. Interface boundaries, 
including the original hadron bubbles as well as any newly formed quark nuggets, invariably evolve from complex 
distorted shapes (influenced by interacting multi-dimensional rarefaction waves, shocks, and phase fronts) to become 
more spherical as the nucleating regions expand, or as the quark nuggets collapse and surface tension becomes more 
important. 

To summarize, although these calculations exhibit complex behavior, there are no signs of hydrodynamic mixing 
instabilities during this transition period, at least for the choice of parameters, scalar field interaction potential, 
equation of state, and grid resolutions investigated. Interfaces between quark and hadron phases remain regular, 
as perturbations along and across the phase boundaries are consistently damped out in time. Our results are thus 
generally consistent with the standard picture of cosmological phase transitions and with the conclusions of Huet 
et al. U, who suggest that electroweak transitions are stable according to linear analysis, even if similar definitive 
statements can not be made of quark-hadron transitions. 

We also note an interesting deflagration 'instability' or acceleration mechanism evident in the third case for which 
we assume an initial thermal discontinuity in space separating different regions of nucleating hadron bubbles. The 
passage of a rarefaction wave through a slowly propagating deflagration can significantly accelerate the condensation 
process and the phase front to velocities near, or in excess of the sound speed. This suggests that if the universe 
super-cools at substantially different rates within causally connected domains, the dominant modes of condensation 
may be through supersonic detonations or fast moving (nearly sonic) deflagrations, assuming that conditions for 
triggering the instability are typical. A similar speculation was made by Kamionkowski and Freese who suggested 
that deflagrations become unstable to perturbations and are converted to detonations over scales determined by 
linear theory and surface tension. In this scenario, instabilities distort the bubble shape, thereby increasing the 
surface area of the wall which accelerates the transition by increasing the rate of condensation and the effective 
velocity of expansion. However, for the calculations presented here, deflagrations are accelerated predominately not 
from turbulent mixing and surface distortion, but from enhanced super-cooling by rarefaction waves generated across 
thermal discontinuities. This effect can be significantly more pronounced in the presence of both longitudinal and 
transverse perturbations as we found by comparing one and two dimensional calculations with the same initial mean 
thermal states. In higher dimensions, the acceleration mechanism is exaggerated further by upwind phase mergers 
due to transverse flow, surface distortion, and mode dissipation effects, a combination that can result in supersonic 
front propagation speeds. 
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